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Vesicle deformations by clusters of transmembrane proteins 
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We carry out a coarse-grained molecular dynamics simulation of phospholipid vesicles with transmembrane 
proteins. We measure the mean and Gaussian curvatures of our protein-embedded vesicles and quantita- 
tively show how protein clusters change the shapes of their host vesicles. The effects of depletion force and 
vesiculation on protein clustering are also investigated. By increasing the protein concentration, clusters are 
fragmented to smaller bundles, which are then redistributed to form more symmetric structures corresponding 
to lower bending energies. Big clusters and highly aspherical vesicles cannot be formed when the fraction of 
protein to lipid molecules is large. 



I. INTRODUCTION 

Biological membranes are found in various com- 
plex shapes^'^ which are closely related to membrane 
functions such as biconcave shape of erythrocytes or 
corkscrew shape of spirochetes. Shape variety depends 
mainly on protein concentration, operation of external 
forces on membranes in cellular environments^, mem- 
brane movement, fusion and budding processes, varia- 
tions in lipid composition, and vesicle trafficking^. Mem- 
brane local curvature is representative for its shape, 
and proteins are believed to play an important role 
in membrane conformation. Proteins behave as both 
generating^ii^ and sensing^ elements for the membrane 
curvature^i^. While protein aggregation induces shape 
transformation, membrane curvature may also generate 
feedback on protein aggregation and yield attractive^i^ 
or repulsive^i^ curvature-mediated interactions between 
them. It has also been shown that bending rigidity and 
membrane thickness affect protein functioning^^^^. 

Proteins affect membrane curvature through various 
ways such as scaffold and local curvature mechanisms, 
and by their integration (as transmembrane proteins) 
with the membrane^i^i^iii. However, the experimental 
measurements of variations in the membrane curvature 
are not easy. The role of proteins on the membrane 
deformations has been studied using continuum elastic 
modeling^i^i^, particle-based^^ and mesoscopio^ simu- 
lations, and hybrid elastic-discrete particle modelsi^. In 
all these studies, lipid bilayers represent a liquid envi- 
ronment with freely diffusing lipid chains that dissolve 
topological deformations. The aggregation of somewhat 
rigid proteins remarkably reduces the diffusion of lipid 
chains and causes shape variations. 

In this study, we use a coarse-grained modeli^ and 
generalize the method of Markvoort et al.^^ to generate 
phospholipid vesicles from initially rectangular and flat 
bilayers with different concentrations of transmembrane 
proteins. We discuss our simulation method in section 
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mi and model the vesicle surface using spherical harmon- 
ics. In section IIIIl we present a method for computing 
the mean and Gaussian curvatures at the locations of 
hydrophilic heads of lipid chains and proteins. The lo- 
cal curvature information together with the sizes of pro- 
tein clusters help us to investigate the deformations of 
host vesicles in section HVl We also study the size distri- 
bution, formation and fragmentation of protein clusters. 
We summarize our fundamental results in section IVl 



II. MODEL DESCRIPTION 

Mesoscopic models have been widely used to study the 
physics of membranes^. Lipids can spontaneously ag- 
gregate and form various membrane^i^U. When lipid 
chains are assembled in the form of a closed 3D surface 
and trap water molecules, a vesicle is generated. En- 
tropy is the main driving mechanism of this proces^. 
To construct vesicles, we insert initially flat rectangu- 
lar lipid bilayers (which may contain proteins) inside a 
box of water molecules. Such a configuration is unstable 
because the tails of boundary lipids are repelled by wa- 
ter molecules and the bilayer is compressed by in-plane 
forces. Consequently, the bilayer buckles and closes itself 
to acquire a minimum potential energy state. Vesicles 
formed through this bilayer vesicle transition process, 
with the progenitor bilayer being surrounded by solvent 
particles (without touching simulation boundaries), have 
a more relaxed pressure distribution. Furthermore, using 
bilayer vesicle transition process we obtain vesicles of 
different sizes in a more controllable process. 

Our simulation box has dimensions of 
with the X-axis being normal to the initial bilayer mid 
plane. We use periodic boundary conditions, and choose 
a sufficiently large box so that the bilayer does not reach 
to the boundaries before vesiculation. The number of 
protein and lipid molecules are denoted by Np and Ni^ 
respectively. The particles that constitute the elements 
of our setups either are water particles (type 1), or have 
hydrophobic (type 2) and hydrophilic (type 3) natures. 
We follow reference [19], and model each lipid chain by 
one hydrophilic head and four hydrophobic tail particles 
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(a) (b) 




FIG. 1. (a) A flexible lipid chain with one hydrophilic head 
(blue) and four hydrophobic tail particles (red), (b) Protein 
molecules consisting of seven strands. Each strand has two 
hydrophilic ends (yellow particles) and five hydrophobic par- 
ticles in the middle (green). Particles of neighboring strands 
are connected by elastic rods (black lines). 



(Fig. [T]a). All particles interact through Lennard- Jones 
(LJ) potential^^ with a cut off radius Vc = 2.5(7. 

We use the integer subscripts i and j for particle 
types, and index a particle in a molecule by the sub- 
script s. For the pairs (i^j) = (1,2) and (2,3) the inter- 
action potential is Vij = 4eij(<Jij/r)^, and for other pairs 
we use Uij = ^eij[[crij /rY'^ - [(Jij/r)^] where crij,eij=l 
{i^j = 1,2,3)^. A harmonic bond potential of the form 
^s^s+i ~ ^6(ks,s+i| — cr)^ is applied between neigh- 
boring particles inside a lipid chain, and we have set 
K}j = 5000es,s+icr~^ where es^s-\-i = 1- With these as- 
sumptions, less than 10% of bond lengthes fluctuate more 
than 2% around a [19]. In some case studies, to include 
bending rigidity in our lipid and protein chains, we use 
the potential 

^s-l,s,s+l — ^bend I J- 7" fl^ \^ 

= i^bend(l - COS(/)s), (2) 

where i^bend is the bending coefficient of non-flexible 
chainsi^. (j)s defines the bending angle between adjacent 
bonds in a single chain. Protein molecules (Fig. [If)) are 
composed of seven strands in a hexagonal arrangement^^. 
A particle in each strand interacts with all of its neigh- 
bors, within the same protein, through the potentials 
U^^,^^ and U^li%^,^^ defined above. In this way, lipid- 
protein interactions are the same as lipid-lipid interac- 
tions. 

In this study, we work with simple proteins whose 
lengths are adjusted in a way that the hydrophobic mis- 
match effect is minimum, and such that a single protein 
does not affect the membrane conformation considerably. 
Shape variations that we report, are thus caused only by 
clustering. Both flexible and rigid chains can be used 
in lipid and protein molecules. The experiments of sec- 
tion HVl show that adding bending rigidity does not in- 
duce remarkable qualitative or quantitative changes in 
the protein tilting angle or in the aggregation products. 
In fact, the strong harmonic bond potentials assumed be- 



tween neighboring strands (in a protein molecule) provide 
enough bending rigidity for our relatively short protein 
molecules. 

We implement equilibrium molecular dynamics simu- 
lation of an NVT ensemble^^ with velocity Verlet algo- 
rithm for integration in the time domain. We use the 
integration time step St = O.OOSto with 

/ TTicr'^ 1 

and set the number density of particle^ to p = 2/3. All 
Particles have equal masses of m = 1, and all lengths are 
scaled by a. The parameters used in our models produce 
correct physical properties of bilayers, including diffu- 
sion coefficients, density profiles and mechanical proper- 
ties like surface tension and stress distributio n^^! Af- 
ter vesicle formation, the position vectors of particles are 
measured with respect to the vesicle center, and the lipid 
or protein heads exposed to water molecules outside the 
vesicle are tagged as surface particles. We assign an in- 
teger number n to each surface particle, and denote the 
total number of surface particles by A^. It is remarked 
that there is not a meaningful correlation between the 
physical location of each particle and its number n. The 
identifier n is used only for statistical purposes. Further- 
more, a single number n is assigned to each surface par- 
ticle, i.e., there are, respectively, one and seven particle 
identifiers corresponding to each lipid chain and protein 
molecule. 



III. LOCAL MEAN AND GAUSSIAN CURVATURES 

To measure the local curvature of vesicles, with and 
without proteins, we express the radial distance of sur- 
face particles from the vesicle center in terms of spherical 
harmonics^l as 

r{e, ^) = E E [«r^r(0, 0) + brBr{<t>, e)] , (4) 

^=0 m=0 

where 

A]^ = Re [Y/^] , BJ^ = Im [Y/^] , (5) 

ym ^ + 1)G " ^) ' (^^^ ^) ^ (g) 

^ ^ ^ y 47r(/ + m)! ^ ^ ^ v / 

Here PJ^ are associated Legendre functions and i = 

We find the coordinates of surface particles (r^, (/>n5 ^n), 

define 

R=[ri{^uOi) r2((/)2,^2) ... r^((/>iv, ^iv) ] ^ ,(7) 

and collect all constants coefficients af^ and b]^ in single 
column vectors a and 6, respectively. The superscript T 
means transpose. We also define the matrices A and B 
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whose elements are 74[^(^^,6>n) and B]^{(t)n-,On)-, respec- 
tively. The discrete form of equation dU thus becomes 



R = Y X, Y=[AB] 



X -- 



(8) 



which in practice, has more equations than unknowns. 

We calculate x using the singular value decomposition 
of Y, and obtain the position of any surface particle from 
dl]) and 

sin(6>)cos((/))^ + sin(6>)sin((/))3*+cos(6>)A^ , (9) 

where {i^j^k) are the unit vectors in Cartesian coordi- 
nates. Let us define and rg as the first order, and r^^, 
r^g*, and rgo as the second order partial derivatives of r 
in ([9j) with respect to (j) and 0. The coefficients of the first 
fundamental form of the surface are thus determined as 



E ■ 



G 



re • re. 



(10) 



, • r^, I' = r^' 

Defining the unit vector normal to the vesicle surface as 
^ = {rcf) X re)/\r(f) x re\, one finds the coefficients of the 
second fundamental form: 

L = r^^- h, M = r^e • N = ree • n. (11) 

The mean curvature H and the Gaussian curvature K at 
the location of each surface particle can thus be computed 
usmg^ 



H 



LN-M^ 
EG-F^ 



K 



EN - 2FM + GL 



(12) 



2{EG-F'^) 

The principal curvatures (Ci,C2) are related to H and 
K through H = {Ci^ C2)/2 and K = CiC'^. In our 
numerical experiments we have truncated the series (|4]) 
at /max = 5. Including /max > 5 terms had ~ 3% im- 
provement in fractional errors. To perform a global shape 
classification of vesicles, we compute the average curva- 
tures 

^ ' H{n) 
K{n) 



(13) 



and their corresponding standard deviations H and K 
for particles living on the surface of model vesicles. 



IV. SIMULATION RESULTS 

We first consider the case of lipid and protein chains 
with no bending rigidity. Our first model, which is called 
1^1, is a vesicle without proteins and it is composed of 
Ni = 2300 similar lipid chains. The initial fiat bilayer 
is bent inside the water (solvent) particles and gradually 
forms an almost spherical vesicle shown in the left panels 
of Fig. [21 To measure the sphericity of this vesicle quan- 
titatively, we use a spherical harmonics expansion with 
/max = 5, and fit a 3D surface to the surface particles. 
The mean local curvature H{n) is computed from ([12]) 
and used in to find H = 0.0564 and H = 0.0068. 
The small value of H compared to H confirms the spher- 
ical nature of vesicle Vi. 



TABLE I. The protein and lipid content, and averaged cur- 
vatures of the simulated vesicles V1-V4. The quantity pb has 
been computed for the bilayers B2-B4. 





Vi 


V2 


V3 


V4 




2300 


2120 


1940 


1760 


Np 





20 


40 


60 


Hp 


_ 


0.0460 


0.0381 


0.0418 


Si 


0.0564 


0.0549 


0.0582 


0.0567 


H 


0.0564 


0.0541 


0.0544 


0.0527 


H 


0.0068 


0.0135 


0.0155 


0.0158 
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0.0030 


0.0029 


0.0029 


0.0027 
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0.0014 


0.0015 


0.0015 


0.0015 
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A. Protein-embedded vesicles 

We replace some lipid chains of our initially fiat bi- 
layer by protein molecules while keeping the area of gen- 
erating bilayer almost constant. We randomly distribute 
proteins in the bilayer sheet, but observe that they form 
small clusters after the bending of membrane and during 
vesicle formation. The vesiculation process takes from 
^ 5 X lO^St for model to 2 x 10^6t for protein- 
embedded vesicles. That is because proteins (or clus- 
ters of proteins) resist against buckling by decreasing the 
fiuidity of the progenitor membrane. After vesicle for- 
mation, we have waited for about 5 x 10^ time steps to 
ensure that vesicles have reached to an equilibrium con- 
dition so that the number and area of protein clusters 
remain constant. 

We have constructed three vesicles with initial bilayers 
of approximately equal surface areas and different pro- 
tein concentrations. We name these vesicles V2, V3 and 
V4, which have been formed inside ^ 161000 water 
molecules. The numbers of lipid and protein molecules 
in our models have been given in Table [H Fig. [2] demon- 
strates three dimensional views of these vesicles and their 
cross sections with the largest diameter so that the ma- 
jority of clusters are visualized. We have also shown some 
water particles inside and outside vesicles. It is seen that 
larger protein clusters have induced lower curvatures in 
their neighborhood, and consequently, prominent defor- 
mations in their host vesicles. 

We also investigate several protein-embedded fiat bi- 
layers where the aggregation of proteins is caused mainly 
by the depletion force. Comparing the sizes and popula- 



FIG. 2. Three dimensional views (top row) and the most informative cross sections (bottom row) of our simulated vesicles 
without (model Vi) and with transmembrane proteins. From left to right: models Vi, V2, V3 and V4. 




FIG. 3. Top views of our flat bilayers B3 (left panel) and B4 
(right panel) with transmembrane proteins. The bilayers B3 
and B4 have the same number of proteins and areas of the 
vesicles V3 and V4, respectively. 



tion of clusters formed in bilayers and vesicles helps us 
to better understand the roles of entropy- and curvature- 
driven aggregation of proteins during vesiculation. We 
simulate three flat bilayers, which extend to the sides of 
the simulation box. The surface areas of these bilayers 
and the population of their proteins match those of the 
vesicles V2-V4. We label these bilayers as ^3 and 
^4. They remain in a flat equilibrium state because of 
periodic boundary conditions that keep them in touch 
with the sides of the simulation box. Fig. [3] displays the 
snapshots of Bs and B4. 

We split the molecules of the outer layer of vesicles 
into two groups of lipids and proteins, and denote by 
Hi{n) and Hp{n) the mean curvatures at the locations 
of lipid and protein heads, respectively. We track only 
initially tagged particles living on the outer surface of 
vesicle (because the probability of flip-flop motions is 
low) and compute the local curvature having their coor- 



dinates (r^, 4>n-,^n) for ^ = 1, 2, . . . , A^. We have plotted 
and {ji = l,p) in Fig. Hfor models F1-V4. The 
curvatures at different surface points have been marked 
by different symbols and colors. The scattered plots pro- 
vide a quantitative insight to the effect of proteins on the 
vesicle shape. It is seen that the fluctuations of H are 
higher in models with proteins. As we mentioned before, 
there is no correlation between the location of particles 
and their corresponding identifier n. Therefore, in Fig. 
[31 points with close values of n are not close physically. 
The oscillatory behavior of the graphs could change by 
renumbering the surface points but the major minima 
always coincide with protein clusters, which flatten their 
host vesicles locally. 

We define Hi and Hp as the averages of mean curva- 
tures (taken over the particles of the same kind) corre- 
sponding to lipid and protein heads, respectively. The 
average mean curvature H for all surface particles (pro- 
tein and lipid heads) and its corresponding standard de- 
viation H have been given in Table |I] together with Hi 
and Hp. Lower values of H for models V2-V4 confirm 
that transmembrane proteins reduce the averaged mean 
curvature through creating low-curvature clusters. Dur- 
ing vesiculation, protein molecules aggregate and form 
clusters. The number of clusters and the population of 
proteins in each cluster are determined by (i) the ran- 
dom arrangement of proteins in the generating bilayer 
membrane (ii) the interplay between depletion force and 
curvature induced interactions (iii) system temperature 
and the concentration of protein molecules. If protein 
molecules are scattered in a vesicle, the membrane will 
maintain its sphericity. However, if the same number of 
proteins build a single cluster, the highest variation in 
vesicle shape will be observed. In our simulations we 
have not seen these two extreme cases. Several clusters 
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FIG. 4. Scattered circles and squares mark, respectively, the mean curvature H (left panels) and the Gaussian curvature K 
(right panels) at the location of proteins (filled red circles) and lipid heads (blue squares) . The averaged curvatures Hp and Hi 
as well as the average Gaussian curvature K and its standard deviation K have been given in Tabled From top to bottom: 
models Vi, V2, V3 and V4. The horizontal axis variable is the normalized real number q — n/N for the nth surface particle. 



with different populations of proteins are usually formed 
in our vesicles. 



B. Curvature-mediated clustering 

Let us denote the number of scattered free proteins of 
a model by Q/, and the number of its clusters by Qc 
Moreover, we indicate by p{m) the number of proteins 
in the mth cluster. For models F2, and V4, we have 
reported the values of Qf, Qc and p{m) in Table HI It is 
seen that the biggest cluster have been formed in vesicle 
V3 and not in V4, which has more proteins. Consequently, 
Hp is lower in vesicle V3 than V4 . The standard deviation 



H is ^ 10% of H for the near-spherical vesicle Vi, but it 
increases remarkably for protein-embedded vesicles. As 
clusters grow, the curvature associated with the ensemble 
of lipids increases, and the vesicle becomes more aspher- 
ical. This can be understood from the larger value of 
Hi in vesicle Vs. A reverse phenomenon is also possible: 
if during the vesicle formation the curvature decreases in 
certain regions, proteins will migrate there and form low- 
curvature clusters. To show the correlation between the 
number of proteins in clusters and H^ we have computed 
the quantity 



Qf^Qc 



Qf 



Qc 

^p(m) 



m—l 



(14) 




TABLE II. The number and sizes of protein clusters for vesi- 
cles V4-V6 with Np = 60 and different initial random arrange- 
ments of proteins. 
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FIG. 5. Three dimensional views of dissociated clusters in 
vesicle Vr with Np = 20 protein molecules initially located in 
a single big cluster. 



C. Convergence tests 



and given its magnitude in Table HI We also define pb 
using a formula similar to ([14]), but for our model bilay- 
ers Bi that correspond to vesicles Vi {i = 2,3,4). The 
computed values of ^5 are given in Table HI 

Since bilayers remain flat during simulation, the ef- 
fect of membrane curvature on the aggregation process 
is ignorable and indicates the contribution of the de- 
pletion force to cluster formation. Comparing the values 
of p and pi) clearly shows that curvature induced inter- 
actions, during vesiculation, facilitate the cluster growth 
and we get p > pb. Moreover, pb is a monotonic function 
of Np while p is not. The reason is the lack of an effective 
fragmentation mechanism in flat bilayers: larger protein 
concentrations always lead to bigger clusters. In vesicles, 
however, the tendency to form a structure with minimum 
bending energy leads to the fragmentation of big clusters 
at the turning (maximum) point of the function p{Np). 

For all surface points including lipid and protein heads, 
we have also calculated (see Table H]) the average Gaus- 
sian curvature K and its standard deviation K using 
(fT2|) and ([l3|). Right panels of Fig. g] show the dis- 
tribution of K for lipid and protein heads. Based on 
the Gauss-Bonnet theorem, any compact manifold A4 
without boundary, is topologically equivalent to a sphere 
and the surface integral K dA of Gaussian curvature 
will be invariant. Since our vesicles become aspherical 
through the shape transformations induced by protein 
clustering, Gauss-Bonnet theorem applies and all mod- 
els V1-V4 must be topologically equivalent. Given that 
the head groups of proteins and lipids are identical in 
our simulations and the areas of generating bilayers are 
almost the same, the surface integral will be approxi- 
mately equal to AK. Data in Table |I] show that both 
K and K remain invariant (from one vesicle to another) 
within a reasonable error threshold. This confirms the 
self-consistency of our models and the results obtained 
from spherical harmonic expansions. 



We have continued our simulations until the size and 
number of clusters become constant in a relaxed equilib- 
rium state. To assure that simulated vesicles are in equi- 
librium, we have carried out various experiments. In the 
first experiment we generated two vesicles from the same 
progenitor bilayer of V4, but using two different sets of 
randomly distributed proteins. We have given the prop- 
erties of new vesicles and Vq in Table [III Although the 
vesicles V4, and Vq start from different initial condi- 
tions, there are minor differences between the properties 
of their clusters. Notably, they posses the same shape 
indicator p=4.61. 

To demonstrate that we usually reach a physical equi- 
librium and not a kinetically trapped state, we designed 
a second experiment and produced a vesicle from an ini- 
tial bilayer where Np=20 proteins (similar to vesicle V2) 
had formed an initial big cluster. In the resulting vesicle 
V7, the proteins of the initial single cluster are dissoci- 
ated into three smaller separate clusters as shown in Fig. 
[5l This result is, again, consistent with the general fea- 
tures of V2 , which had been obtained from a completely 
different initial condition. It is worth noting that we 
have observed a transient interplay between clustering 
and fragmentation well before reaching the equilibrium 
state. 

We have repeated our simulations with p = 0.8 and 
with non- flexible lipid and protein chains. Using a bend- 
ing stiffness i^bend = 5 in lipids and i^bend = 80 in pro- 
tein strands, the bilayer vesicle transition is slowed 
down but we do not observe considerable change, either 
quantitative or qualitative, in the clustering phenomenon 
and the shape transformation of vesicles: the numbers 
and sizes of final clusters and the shape parameter p are 
similar in all models. By making stiffer molecular chains 
and increasing the density, lipid diffusion is decreased, 
which in turn, yields a longer relaxation time. 
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V. CONCLUSIONS 

We have studied the phenomena of protein clustering 
and membrane shape transformation during membrane 
vesiculation and afterwards. Comparing relaxed vesicles 
and bilayers shows that protein clustering during vesicu- 
lation occurs due to both entropy-driven depletion force 
and the curvature- mediated interactions. The latter ef- 
fect enhances the generation of larger protein clusters 
and determines bilayer's bending rigidity. Once the vesi- 
cle is formed, protein clusters locally flatten their host 
vesicles and increase the bending energy as H and Hi 
increase. The system, however, cannot tolerate the in- 
crease in the bending energy for protein concentrations 
beyond a critical value. By increasing the protein con- 
centration, bigger protein clusters are not formed in our 
simulations, or they break apart. We anticipate a uni- 
form distribution of fragmented clusters, like the shape 
of a soccer ball. Our observations show that low protein 
concentrations do not lead to efficient cluster formation. 
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